% !TEX root = FDS_Technical_Reference_Guide.tex

% \usepackage{tikz,tikz-3dplot}
% \usetikzlibrary{arrows}

% \newenvironment{myfont}{\fontfamily{\ttdefault}\selectfont}{\par}

% \mathchardef\mhyphen="2D


\typeout{new file: Complex_Geometry_Chapter.tex}


\chapter{Unstructured Geometry (Beta)}
\label{sec:unstructured_geometry}

Unstructured geometry is treated using a ``cut-cell'' (CC) method, following the seminal work of Michael Aftosmis and Marsha Berger \cite{Berger:2012,Berger:2017,May:2017}.  The term cut-cell refers to the gas phase region of a Cartesian grid cell that has been ``cut'' by the unstructured geometry partitioning the cell into gas and solid sub-volumes.  Faces of the cut-cell that have also been carved from Cartesian faces and connected to other cut-cells are called \emph{gasphase cut-faces} (these are gas-gas interfaces).  The faces on the boundary of the solid in a cut-cell are \emph{in-boundary cut-faces} (these are gas-solid interfaces).  There are also faces of a cut-cell which are not cut and which connect the cell to a neighboring regular Cartesian gas phase cell; these faces are referred to as \emph{regular faces} (also a gas-gas interface).

The introduction of cut-cells surrounding an embedded boundary generates the need for an unstructured solution algorithm.  In particular, the pressure equation is unstructured and so fast trigonometric solvers are not viable (i.e. the \ct{FFT} Poisson solver).  Instead, an unstructured local matrix (\ct{ULMAT}) solver has been developed using the Pardiso and Sparse Cluster Solvers available in the Intel Math Kernel Library (MKL).

The following sections describe the finite volume methods used to discretize the governing equations in the unstructured cut-cells surrounding embedded boundaries.

\section{Scalar Transport Discretization for Complex geometry}

\subsection{Multispecies Mass Balance Equations}

Consider a set of gaseous, reacting  chemical species $\alpha=1,\dots,N$ flowing on a given spatial domain $\Omega \in \mathbb{R}^n, \; n=2,3$, with boundary $\partial \Omega$, paramaterized by an Eulerian reference frame $N$. These species are transported on a given point $\mathbf{x}$ in space with velocity $\mathbf{u}_\alpha(\mathbf{x},t)$ respect to $N$, and a mass weighted average velocity $\mathbf{u}(\mathbf{x},t)$
%
\begin{equation}
  \mathbf{u} = \frac{ \sum\limits_{\alpha=1}^{N} {\rho_\alpha \mathbf{u}_\alpha}}{\rho} \; , \; \rho =  \sum\limits_{\alpha=1}^{N} {\rho_\alpha} \label{eq:veldens}
\end{equation}
%
where space and time dependencies are not shown for simplicity, $\rho_\alpha(\mathbf{x},t) = \rho(\mathbf{x},t) Y_\alpha (\mathbf{x},t)$, $\rho$ is the mixture density and $Y_\alpha = \rho_\alpha / \rho$ is species $\alpha$ mass fraction. Definition: $\rho_\alpha(\mathbf{x},t)$ is the amount of mass of species $\alpha$ on a given differential volume $d\Omega$ centered in $\mathbf{x}$, while  $\rho(\mathbf{x},t)$ is the mass of \textit{all} species in $d\Omega$.

As defined, each species $\alpha$ has a velocity $\mathbf{u}_\alpha(\mathbf{x},t)$ which is different from the average mixture velocity $\mathbf{u}(\mathbf{x},t)$. This difference quantity is called the \textit{diffusion velocity} of species $\alpha$, $\mathbf{V}_\alpha$
%
\begin{equation}
   \mathbf{V}_\alpha(\mathbf{x},t) = \mathbf{u}_\alpha(\mathbf{x},t) - \mathbf{u}(\mathbf{x},t) \label{eq:vdiff}
\end{equation}
%

The transport-reaction mass balance equation for each individual species $\alpha$ is given by
%
\begin{equation}
   \frac{\partial \rho_\alpha}{ \partial t} + \nabla \cdot (\rho_\alpha  \mathbf{u}_\alpha) = \dot{m}_\alpha''' \; , \; \alpha=1,\dots,N \label{eq:bal}
\end{equation}
%
in $\Omega$, where $\dot{m}_\alpha'''(\mathbf{x},t)$ is the chemical reaction volume source (area source for $n=2$, sink if negative) of species $\alpha$. The problem is completely defined stating initial and boundary conditions for $\rho_\alpha(\mathbf{x},t)$ as well as time dependent fields $\mathbf{u}_\alpha(\mathbf{x},t)$, $\dot{m}_\alpha'''(\mathbf{x},t)$.

Noting that $ \rho_\alpha = \rho Y_\alpha$, and using equation~\eqref{eq:vdiff} in the previous
%
\begin{equation}
   \frac{\partial \rho Y_\alpha}{ \partial t} + \nabla \cdot \left( \rho Y_\alpha  (\mathbf{u}+\mathbf{V}_\alpha) \right) = \dot{m}_\alpha''' \; , \; \alpha=1,\dots,N \label{eq:bal2}
\end{equation}
%

We see from these equations that the evolution of these species can be considered as a function of the mixture $\rho$ and average $\mathbf{u}$, and written as such, depending on the other species through Eqs.~\eqref{eq:veldens}, and through the mass conservation statement on reactive sources
%
\begin{equation}
  \sum\limits_{\alpha=1}^{N} \dot{m}_\alpha'''(\mathbf{x},t) = 0
\end{equation}
%

In Eq.~\eqref{eq:bal2}, the following convective and diffusive fluxes per unit time and area are found:
%
\begin{eqnarray}
  \mathbf{J_{c \alpha}} &=& \rho Y_\alpha  \mathbf{u} \label{eq:jc} \\
  \mathbf{J_{d \alpha}} &=& \rho Y_\alpha  \mathbf{V}_\alpha \label{eq:jd}
\end{eqnarray}
%
The definition of diffusion velocities $\mathbf{V}_\alpha$ is a central component of multispecies mass transport. For atmospheric combustion we consider the binary diffusion of all species respect to Nitrogen, the most abundant (\textit{background}) species. This simplification allows us to uncouple the diffusive fluxes, using Ficks Law of diffusion
%
\begin{equation}
   \mathbf{J_{d \alpha}} = \rho Y_\alpha  \mathbf{V}_\alpha = - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha
\end{equation}
%
where $D_\alpha$ is the binary diffusivity coefficient respect to the background species.

Using this last expression in Eq.~\eqref{eq:bal2} we arrive to the basic form of the species balance equations based on mass fractions, and used in FDS
%
\begin{equation}
   \frac{\partial \rho Y_\alpha}{ \partial t} + \nabla \cdot ( \rho Y_\alpha  \mathbf{u} ) = \nabla \cdot ( \rho D_\alpha \boldsymbol{\nabla} Y_\alpha ) + \dot{m}_\alpha''' \; , \; \alpha=1,\dots,N \label{eq:bal3}
\end{equation}
%

As will be seen in next section, a spatially discretized version of the previous equations in a domain divided in $ii=1,\ldots,n_{tot}$ computational cells is given by:

%
\begin{equation}
\frac{\partial \left(\rho Y_\alpha \right)_{ii} }{\partial t} V_{ii}  = F^{adv}_{\alpha,ii} + F^{diff}_{\alpha,ii} + \dot{m}_{\alpha,ii}''' +  F^{BC}_{\alpha,ii} \; , \; \alpha=1,\dots,N \label{eq:discbal3}
\end{equation}
%
where $V_{ii}$ is the cell volume, $\left(\rho Y_\alpha \right)_{ii}$ are the cell unknown partial densities, 
$F^{adv}_{\alpha,ii}$, $F^{diff}_{\alpha,ii}$, $\dot{m}_{\alpha,ii}''$, $F^{BC}_{\alpha,ii}$ are terms due to advection, diffusion, reaction and boundary conditions. Given that the time integration method we use is explicit, these vectors make use of available values of $\left(\rho Y_\alpha \right)_{ii}$ from a previous time step or substep. 


\subsection{Finite Volume Method on Staggered Grid in 2D}

The finite volume discretization (FV) starts by considering the integral form of Eq.~\eqref{eq:bal3} over a cell control volume $\Omega_{ii}$. For a cut-cell
$ii$ individualized by the index pair $(i,j)$ we have
%
\begin{equation}
 \int_{\Omega_{ii}} {\frac{\partial \rho Y_\alpha}{\partial t}} d \Omega + \int_{\Omega_{ii}} { \boldsymbol{\nabla} \cdot  \left(  \rho Y_\alpha \mathbf{u} \right)
      } d \Omega  = -\int_{\Omega_{ii}} { \boldsymbol{\nabla} \cdot \left(  \mathbf{J_{d \alpha}}  \right)  } d \Omega + \int_{\Omega_{ii}} { \dot{m}_\alpha''' } d \Omega \label{eq:intconvdiff}
\end{equation}
%
Assuming a time independent control volume, the time derivative and source terms are approximated by
%
\begin{eqnarray}
  \int_{\Omega_{ii}} {\frac{\partial \rho Y_\alpha}{\partial t}} d \Omega & = & \frac{\partial}{\partial t} \int_{\Omega_{ii}} {\rho Y_\alpha} d \Omega
  = \frac{\partial \overline{\rho \: Y_\alpha }_{i,j}}{\partial t} V_{i,j} \\
  \int_{\Omega_{ii}} { \dot{m}_\alpha''' } d \Omega & = & \overline{ \dot{m}_\alpha''' }_{i,j} V_{i,j} \label{eq:intcons}
\end{eqnarray}
%
where $V_{i,j}$ is the volume of cell $(i,j)$. These cell averaged quantities (denoted by tildes) match the cell centroid values of the corresponding scalar fields up to second order spatial accuracy. A consequence of this known fact is that using a second order finite difference method (FD), or a \textit{difference} finite volume approach with same definition
 for diffusion flux spatial derivatives and interpolation on the advective term, will lead to identical discretization matrices (they will vary up to a volume factor) on uniform Cartesian grids. The difference between these two numerical discretization methods arises from the source $\dot{m}_\alpha''' $, which for a given cell is the cell centroid value $\dot{m}'''_{\alpha (i,j)}$ in FD, and it is the cell average $ \widetilde{ \dot{m}_\alpha''' }_{i,j}$ in FV. To discretize both diffusive and advective terms in equation~\eqref{eq:intconvdiff}, we make use of the divergence theorem. In the following we drop the tildes to simplify the notation, keeping in mind that using FV, quantities will always be cell or face averaged where it corresponds.


\subsection{Finite-Volume Discretization in Cut-Cells} \label{sec:cc}

The treatment of irregular gas phase Eulerian grid faces and cells that remain after the inclusion of immersed bodies is described in this section.
We assume (see Fig.~\ref{Fig:FVdiscCC}a) the cut-cell definition algorithm has been successful in defining the following:
%
\begin{itemize}
   \item All \texttt{GASPHASE} faces in the $x$ and $y$ directions, which are tagged as regular or cut-faces.
           For these, we also know vertex points, areas (length in 2D), and face centroid location. We know for faces that are \texttt{GASPHASE}
           cut-faces the $i,j$ coordinates of the face they belong to.
  \item All \texttt{INBOUNDARY} faces, which arise from the intersection of the immersed bodies surface elements and the Eulerian
           grid cells. For these we know vertex points, areas (length in 2D), and face centroid location.
   \item All \texttt{GASPHASE} cut cells and regular cells. We know which faces define their boundary. For cut-cells, a list of cut-faces
           (type \texttt{GASPHASE} and \texttt{INBOUNDARY}) is provided. Also we know their volume and face centroid location.
\end{itemize}
%
\begin{figure}[h]
      %\centering
      \includegraphics[trim = 65mm 50mm 70mm 40mm, clip,
       width=0.49\linewidth]{../../../fig/fds/Geometry_Figures/CutCellsSketch.png}
      \includegraphics[trim = 65mm 40mm 50mm 30mm, clip,
       width=0.49\linewidth]{../../../fig/fds/Geometry_Figures/CutCellsBCSketch.png}
      \put(-350,-10){(a)}
      \put(-120,-10){(b)}
      \caption[Cut-cell FV discretization in 2D]{Cut-cell FV discretization in 2D: (a) \texttt{GASPHASE} Cut-cell  $ii$ is bounded by: $\mathbf{1,2}$ \texttt{INBOUNDARY} faces, $\mathbf{3,6}$ \texttt{GASPHASE} cut-faces, and $\mathbf{4,5}$  \texttt{GASPHASE} regular faces.  (b) Geometrical elements used in the discretization of the diffusive term for \texttt{GASPHASE} cut-face $\mathbf{3}$, and \texttt{INBOUNDARY} face $\mathbf{1}$.}
   \label{Fig:FVdiscCC}
\end{figure}
%
\subsection{Cut-Cell Discretization of Diffusive Term}  \label{sec:CCdiff}

Consider the FV discretization of the diffusive term (equation~\eqref{eq:intconvdiff}) on cut-cell $ii$ in Fig.~\ref{Fig:FVdiscCC}b. We have
%
\begin{equation}
    \int_{\Omega_{ii}} { \boldsymbol{\nabla} \cdot \left(  \mathbf{J_{d \alpha}}  \right)  } d \Omega =
    \int_{\partial \Omega_{ii}} { \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right) \cdot \hat{\mathbf{n}}_{ii} } \: d \partial \Omega = \sum^{nf_c=6}_{k=1}
    \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right)_k \cdot \hat{\mathbf{n}}_{ii,k} \: A_k \label{eq:discfvdiffcc}
\end{equation}
%
The $nf_c=6$ faces that compose the boundary of the cut-cell can be divided in:

\subsection*{A. Faces $k=\mathbf{4},\mathbf{5}$ are \textit{regular} \texttt{GASPHASE} Faces Connecting a Cut-Cell with a Regular Cell:}
The treatment of these is as done in regular Cartesian faces, with the caveat that the diffusive flux computation now involves the cut-cell centroid location (i.e. to compute spatial derivatives of $Y_\alpha$ and interpolation of $\rho D_\alpha$ in equation~\eqref{eq:discfvdiffcc}).

\subsection*{B. Faces $k=\mathbf{3},\mathbf{6}$ are \texttt{GASPHASE} \textit{Cut-Faces}:}

In order to compute the discrete term of Eq.~\eqref{eq:discfvdiffcc} on these faces, the factor $\rho D_\alpha$ needs to be interpolated from cell centroids to the face centroid and the spatial derivative of $Y_\alpha$ computed at the face centroid. The use of face centroids is required to maintain spatial accuracy when going from the left to right hand side of Eq.~\eqref{eq:discfvdiffcc} in difference FV methods. To maintain accuracy of the overall cut cell method, interpolation and differentiation to the cut-face centroid must also be sufficiently accurate.

The following discussion is fairly general, in the sense that it is agnostic to the interpolation methods employed (Lagrange polynomials, isoparametric, least squares, etc.).
Consider a stencil of points $e=1,...,ne$, where the scalar $Y_\alpha$ is assumed defined (i.e. cell centroids). Then for a location of interest $ck$ (i.e. the centroid of \texttt{GASPHASE} cut-face $k=3,6$), interpolation and derivatives of $Y_\alpha$ at $ck$ can be obtained as
%
\begin{eqnarray}
   Y_\alpha(\mathbf{x}_{ck}) & = & \sum^{ne}_{e=1} \phi_e(\mathbf{x}_{ck}-\mathbf{x}_e) Y_{\alpha e} \label{eq:interpqfluid} \\
   \frac{\partial Y_\alpha(\mathbf{x}_{ck})}{\partial x_i} & = & \sum^{ne}_{e=1} \frac{\partial \phi_e(\mathbf{x}_{ck}-\mathbf{x}_e)}{\partial x_i} Y_{\alpha e} \label{eq:interpdqfluid}
\end{eqnarray}
%
where $\mathbf{x}_e$ is the location of stencil point $e$, $Y_{\alpha e} =Y_\alpha(\mathbf{x}_e)$  and $\phi_e(\mathbf{x}-\mathbf{x}_e)$, $e=1,...,ne$ is a suitable set of interpolation functions.
For this discussion we assume that no boundary values of $Y_\alpha$ are involved on the interpolation. For these faces the diffusive flux is:
%
\begin{eqnarray}
  \left( - \rho D_q \boldsymbol{\nabla} Y_\alpha \right)_k & = & - \sum^{ne}_{m=1} \phi_m(\mathbf{x}_{ck}-\mathbf{x}_m) [\rho D_\alpha]_m \times
      \sum^{ne}_{e=1} \boldsymbol{\nabla} \phi_e(\mathbf{x}_{ck}-\mathbf{x}_e) \; Y_{\alpha e} \nonumber \\
      &=& - [\rho D_\alpha]_k \sum^{ne}_{e=1}  \boldsymbol{\nabla} \phi_e(\mathbf{x}_{ck}-\mathbf{x}_e) \; Y_{\alpha e}
      \label{eq:diffflx}
\end{eqnarray}
%
where
%
\begin{equation}
 [\rho D_\alpha]_k = \left[ \sum^{ne}_{m=1} \phi_m(\mathbf{x}_{ck}-\mathbf{x}_m) [\rho D_\alpha]_m \right]
\end{equation}
%

%%
\subsection*{1D Linear Isoparametric Polynomial Interpolation:}
%%%%%

Consider a stencil with cells $jj$ and $ii$ as in Fig.~\ref{Fig:FVdiscCC}b. Their centroids are located at $\mathbf{x}_{jj}$ and $\mathbf{x}_{ii}$ respectively. Parameterizing point locations along the segment that unites $\mathbf{x}_{jj}$ to $\mathbf{x}_{ii}$ by $-1 \le \xi \le 1$, we have
%
\begin{eqnarray}
   \mathbf{x}(\xi)  & = &  \phi_1(\xi) \mathbf{x}_{jj} + \phi_2(\xi) \mathbf{x}_{ii} \label{eq:xxi} \\
   Y_\alpha (\xi)  & = & \phi_1(\xi) Y_{\alpha jj} + \phi_2(\xi) Y_{\alpha ii}
\end{eqnarray}
%
\begin{equation}
[ \rho D_\alpha ] (\xi)  =  \phi_1(\xi) [ \rho D_\alpha ]_{jj} + \phi_2(\xi) [ \rho D_\alpha ]_{ii}
\end{equation}
where
%
\begin{equation}
   \phi_1(\xi) = \frac{1}{2} \left(1 - \xi \right) \; ; \; \phi_2(\xi) = \frac{1}{2} \left(1 + \xi \right)
\end{equation}
%
We assume the value of the parameter where the cut-face lies known $\xi_p$, corresponding to the face and centroid segment intersection $\mathbf{x}_p$. From Eq.~\eqref{eq:xxi} with $\mathbf{x}=x\hat{\mathbf{i}}+y\hat{\mathbf{j}}$, we have
%
\begin{eqnarray}
   \frac{\partial \xi}{\partial x} = \frac{2}{(x_{ii}-x_{jj})} \; ; \;  \frac{\partial \xi}{\partial y} = \frac{2}{(y_{ii}-y_{jj})}
\end{eqnarray}
%
and the interpolation function gradient components are:
%
\begin{eqnarray}
  \frac{\partial\phi_1}{\partial x} & = & \frac{\partial\phi_1}{\partial \xi} \frac{\partial \xi}{\partial x}=-\frac{1}{(x_{ii}-x_{jj})} \\
  \frac{\partial\phi_1}{\partial y} & = & \frac{\partial\phi_1}{\partial \xi} \frac{\partial \xi}{\partial y}=-\frac{1}{(y_{ii}-y_{jj})} \\
  \frac{\partial\phi_2}{\partial x} & = & \frac{\partial\phi_2}{\partial \xi} \frac{\partial \xi}{\partial x}= \frac{1}{(x_{ii}-x_{jj})} \\
  \frac{\partial\phi_2}{\partial y} & = & \frac{\partial\phi_2}{\partial \xi} \frac{\partial \xi}{\partial y}= \frac{1}{(y_{ii}-y_{jj})}
\end{eqnarray}
%
which are constant. Finally, being $k=3$, for cut-cell $jj$ the normal is $\hat{\mathbf{n}}_{jj,k}=\hat{\mathbf{j}}$, and expression~\eqref{eq:diffflx} reduces to
%
\begin{equation}
  \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right)_k \cdot \hat{\mathbf{n}}_{jj,k} \: A_k =
  -[ \rho D_\alpha ]_p \frac{\left(  Y_{\alpha ii} - Y_{\alpha jj} \right)}{(y_{ii}-y_{jj})}A_k
\end{equation}
%
Note that, an assumption of approximation of the cut-face centroid location by point $p$ has been made. For cut-cell $ii$ on the high side of face $k=3$ the expression is
%
\begin{equation}
  \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right)_k \cdot \hat{\mathbf{n}}_{ii,k} \: A_k =
  [ \rho D_\alpha ]_p \frac{\left(  Y_{\alpha ii} - Y_{\alpha jj} \right)}{(y_{ii}-y_{jj})}A_k
\end{equation}
%

Similar expressions are found for a \texttt{GASPHASE} cut-face normal to the $x$ direction. 


\subsection*{C. Faces $\mathbf{1},\mathbf{2}$ are \texttt{INBOUNDARY} \textit{Cut-Faces}:}

% Here define how we treat CFACEs Boundary Conditions.
Consider \texttt{INBOUNDARY} cut-face $k=\mathbf{1}$ in Fig.~\ref{Fig:FVdiscCC}b. If a Neumann boundary condition is specified in such a face, the term
%
\begin{equation}
  \left( - \rho D_\alpha \boldsymbol{\nabla} Y_\alpha \right)_k \cdot \hat{\mathbf{n}}_{ii,k} \: A_k = [\rho D_\alpha]_k \partial \gamma_{ii,k} \: A_k
\end{equation}
%
where as before, the prescribed normal derivative of $Y_\alpha$ on this face $\partial Y_\alpha / \partial n |_k=\partial  \gamma_{ii,k}$ is defined in the direction of $\hat{\mathbf{n}}_1$.
This term adds to $\mathbf{F_{BC \alpha}}$ in the location of cell $ii$. For \texttt{SOLID\_BOUNDARIES}, zero diffusive mass flux implies a homogeneous $\partial Y_\alpha / \partial n |_k=0$ for the mass fraction field of species $\alpha$.
In case a Dirichlet boundary condition is specified in boundary face $k=1$ the normal mass fraction derivative is taken as:
%
\begin{equation}
   \frac{\partial Y_\alpha}{\partial n} |_{k=1}=\frac{Y_{Ec1} - Y_{\alpha c1}}{h} 
\end{equation}
%
where the value of $Y_{Ec1}$ is the cell value, and $h$ is taken as half the cell volume divided by surface boundary area on the cell. Note that, as discussed later in Sec.~\ref{sec:exscl}, when small cells are linked to neighboring cells these quantities refer to the control volume of linked cells.


\subsection{Cut-Cell Discretization of Advective Term}

We study the discretization of the advective term for the six faces of cell $ii$ in Fig.~\ref{Fig:FVdiscCC}b:
%
\begin{equation}
  \int_{\Omega_{ii}} { \boldsymbol{\nabla} \cdot  \left(  \rho Y_\alpha \mathbf{u} \right) } d \Omega =
  \int_{\partial \Omega_{ii}} { \left( \rho Y_\alpha \mathbf{u} \right) \cdot \hat{\mathbf{n}}_{ii} } \: d \partial \Omega =
  \sum^{nf_c}_{k=1} \left(  \overline{[\rho Y_\alpha]} \mathbf{u} \right)_k \cdot \hat{\mathbf{n}}_{ii,k} \: A_k \label{eq:fvadv}
\end{equation}
%
where, as required for stability in the treatment of the advective term, $\overline{[\rho Y_\alpha]}_k$ represents a flux limited interpolation to the face. See Sec.~\ref{sec_flux_limiters}. For completeness, the use of interpolation functions in the interpolation process to a face, can also be extended to flux limited interpolation. Consider for a generic face $k$ the advective flux for the conserved quantity $\rho Y_\alpha$ off cell $ii$:
%
\begin{equation}   
    \left(  \overline{[\rho Y_\alpha]} \mathbf{u} \right)_k \cdot \hat{\mathbf{n}}_{ii,k}  = \overline{[\rho Y_\alpha]}_k \left(  \mathbf{u}_k \cdot \hat{\mathbf{n}}_{ii,k} \right) 
\end{equation}
%
Then, a flux limited spatial interpolation of the form~\eqref{eq:interpqfluid} for $[\rho Y_\alpha]$ to face centroid $ck$ is
%
\begin{equation}
  \overline{[\rho Y_\alpha]}_k = \sum^{ne}_{e=1} \overline{\phi}_e(\mathbf{x}_{ck}-\mathbf{x}_e, \mathbf{u}_k) \left( \rho Y_{\alpha} \right)_e
\end{equation}
%
where $e=1,...,ne$ refers to a suitable flux limited stencil of \texttt{GASPHASE} cell centroids. The $\overline{\phi}_e$ are the flux limited interpolation functions  related to said centroids.

\subsection*{A. Faces $k=\mathbf{4},\mathbf{5}$ are \textit{regular} \texttt{GASPHASE} Faces Connecting a Cut-Cell with a Regular Cell:}
The treatment of these is as described in the previous section, with the caveat that the advective flux computation for $\rho Y_\alpha$ may now involve a four point stencil ($e=1,...,4$) spatial flux matched interpolation to the face. Consider the situations depicted in Fig.~\ref{Fig:FVFlxLimCC} for a regular face $RC$ connecting regular cell and cut-cell. In the left figure both extreme interpolation points of the stencil have associated species densities $\rho Y_{\alpha 1}$ and $\rho Y_{\alpha 4}$, as at least one cut-cell or regular cells are defined at those locations. In the right figure the lower extreme position of the stencil (cell S) lies completely within the \texttt{SOLID} phase, and it takes the value of the inner cell $\rho Y_{\alpha 2}$, rendering the interpolation to RC of Godunov type. Note also the loss of accuracy on the flux limited interpolation described in Sec.~\ref{sec_flux_limiters} due to cut-cell centroids making the stencil not uniformly spaced.
 
\begin{figure}[h]
      %\centering
      \includegraphics[trim = 5mm 75mm 15mm 60mm, clip,
       width=0.95\linewidth]{../../../fig/fds/Geometry_Figures/CutCellsFLXLIMSketch.png}
      \put(-330,-10){(a)}
      \put(-110,-10){(b)}
      \caption[Flux limited interpolation]{Flux limited interpolation: (a)  Regular and cut-cell species densities are available for 4 point stencil flux limited interpolation to face $RC$.  (b) Lower cell in stencil is within solid region (cell S), therefore inner cell value $\rho Y_2$ is used for this extreme point.}
   \label{Fig:FVFlxLimCC}
\end{figure}

\subsection*{B. Faces $k=\mathbf{3},\mathbf{6}$ are \texttt{GASPHASE} \textit{Cut-Faces}:}

For \texttt{GASPHASE} cut-face $k=3$ in Fig.~\ref{Fig:FVdiscCC}b, the advective term corresponding to cell $jj$ is
%
\begin{equation}
  \left( \rho Y_\alpha \mathbf{u} \right)_k \cdot \hat{\mathbf{n}}_{jj,k} \: A_k = \overline{[\rho Y_\alpha]}_k \left( \mathbf{u}_k \cdot \hat{\mathbf{n}}_{jj,k} \right) \: A_k \label{eq:convgcutface}
\end{equation}
%
where the average face normal velocity $\mathbf{u}_k \cdot \hat{\mathbf{n}}_{jj,k} $ (i.e. approximated by the face centroid value) is assumed known. The quantity $\overline{[\rho Y_\alpha]}_k$ is obtained by flux limited interpolation to the face in similar manner as described previously. The flux component for cut-cell $ii$ is naturally opposite in sign to what is computed for cell $jj$.


\subsection*{C. Faces $\mathbf{1},\mathbf{2}$ are \texttt{INBOUNDARY} \textit{Cut-Faces}:}

The treatment of advective terms in boundary cut-faces is similar as described before. Consider a normal velocity $\mathbf{u}_1=v_n \hat{\mathbf{n}}_{ii,1}$ (into solid region) imposed at the boundary in \texttt{INBOUNDARY} cut-face k = 1 in Fig.~\ref{Fig:FVdiscCC}b. The advective component is in this case:
%
\begin{equation}
   \left( \rho Y_\alpha \mathbf{u} \right)_1 \cdot \hat{\mathbf{n}}_{ii,1} \: A_k = \overline{[\rho Y_\alpha]}_1 v_n \: A_k
\end{equation}
and the value of $\overline{[\rho Y_\alpha]}_1=\rho Y_{\alpha ii}$ if $v_n > 0$, or $\overline{[\rho Y_\alpha]}_1=\rho \gamma_{ii,1}$ prescribed otherwise.


\subsection{Unsteady Evolution: Explicit Time Integration for Scalars} \label{sec:exscl}

It is well known that cut-cell methods pose a significant time constraint when used with explicit time integration methods.
In general, there will arise \texttt{GASPHASE} cut-cells whose small-size will severely penalize
the time step. We recall that each cell on the gas phase, including cut-cells, needs to meet CFL and Von Neuman stability constraints. Several different ways have been proposed in the literature to deal with this problem, i.e., cell merging, mixing or linking methods. In general, these lead to ad hoc selection procedures for surrounding cells, having to deal with many special cases, and in some cases potential solution deterioration close to the boundary.

We use a simple and robust procedure to address this problem. Within the scheme for numeration of scalar cell unknowns, a test is performed on cut-cells. If the cut-cell volume is less than the threshold volume $V_{thr}= $\ct{CCVOL\_LINK} $V_{cart}$, where $V_{cart}$ is the local Cartesian cell volume and \ct{CCVOL\_LINK}$<1$ is a threshold factor (default value 0.5), the unknown number this cell takes is the one of an adjacent cell which has a volume larger than $V_{thr}$. This mathematically defines a single control volume of the two linked cells. Cell volumes are added in building the mass matrix for the FV discretization, and fluxes and matrix terms are added with their corresponding signs. Note that, flux quantities corresponding to the common face of these two cells effectively cancel on the single equation for the set.
Alternatively, if after a number of cell numbering iterations (default 2) an unlinked small cell persists in the mesh, this small cell is blocked and its corresponding Cartesian cell tagged as solid. % Default 2 link numbering iterations.

\section{Momentum Time Marching and Immersed Boundaries}

\subsection{Momentum-Pressure Coupling}

As a first approximation, consider the Newtonian flow problem defined by the following set of partial differential equations:
%
\begin{eqnarray}
  \frac{\partial \mathbf{u}(\mathbf{x},t)}{\partial t} &=& - \left[ \mathbf{F}(\mathbf{u},\mathbf{x},t) + \boldsymbol{\nabla} H(\mathbf{x},t) \right] \; , \; \mathbf{x} \in \Omega - \sum{\Omega_i} \; , \; t \in \mathbb{R}_+ \label{eq:LowMachMom} \\
         \nabla \cdot \mathbf{u} (\mathbf{x},t) & = & \left(\nabla \cdot \mathbf{u} \right)^{th} \label{eq:LowMachDiv}
\end{eqnarray}
%
where equation~\eqref{eq:LowMachMom} is the momentum equation, $\mathbf{u}(\mathbf{x},t)$ is the spatial velocity field, $\mathbf{F}(\mathbf{u},\mathbf{x},t)$ is a vector containing convective, diffusive and possibly other force terms, and $H(\mathbf{x},t)$ is a potential scalar field (physically the head field in this case, commonly called pressure). For sake of argument here, to represent the low Mach number approximation employed in FDS, it suffices to consider a specified divergence field $\left(\nabla \cdot \mathbf{u} \right)^{th} (\mathbf{x},t)$ (thermodynamic divergence). This divergence field in the thermally buoyant flow model used in FDS is a proxy for the energy equation.
The domain $\Omega - \sum{\Omega_i}$ represents the fluid region, and boundary conditions are prescribed for $\mathbf{u}(\mathbf{x},t)$ on $\partial \Omega,\partial \Omega_1,...,\partial \Omega_{nbods}$.

Classical fractional step methods for time integration of incompressible or low Mach number flow are based on two operations: First, momentum transport to obtain intermediate velocities, and second, projection of velocities into a target divergence field. Consider the Forward Euler update of the governing equations from $t_n$ to $t_{n+1}=t_n + \Delta t$ of the form: Given $ \mathbf{u}^n=\mathbf{u}(\mathbf{x},t_n)$, $\nabla \cdot \mathbf{u}^{n+1} = \left( \nabla \cdot \mathbf{u}^{n+1} \right)^{th}$ known
%
\begin{eqnarray}
  \frac{\mathbf{u}^{n+1}-\mathbf{u}^{n}}{\Delta t} &=& - \left[ \mathbf{F}^n +  \boldsymbol{\nabla} H^n \right] \label{eq:LowMachMomEu}\\
  \nabla \cdot \mathbf{u}^{n+1} &=& \left( \nabla \cdot \mathbf{u}^{n+1} \right)^{th} \label{eq:LowMachDivEu}
\end{eqnarray}
%
where as time has been discretized, $\mathbf{u}^{n+1}$ represents a numerical approximation to the solution in Eqs.~\eqref{eq:LowMachMom}-\eqref{eq:LowMachDiv} at time $t_{n+1}$. As the potential field $H(\mathbf{x},t)$ does not have a time evolution equation, it is assumed responsible of enforcing the divergence condition and used on the projection step. Taking the divergence of equation~\eqref{eq:LowMachMomEu} and considering  the constraint Eq.~\eqref{eq:LowMachDivEu}, the two steps of the method are
%
\begin{enumerate}
  \item Solve Poisson equation for $H^n$:

\begin{equation}
   \nabla \cdot \boldsymbol{\nabla} H^n = - \left[ \frac{\left( \nabla \cdot \mathbf{u}^{n+1} \right)^{th} - \nabla \cdot \mathbf{u}^{n}}{\Delta t} \right] - \nabla \cdot \mathbf{F}^n \label{it:FSPoisson}
\end{equation}

  \item Obtain final velocity for step:

  \begin{equation}
     \mathbf{u}^{n+1} = \mathbf{u}^{n} - \Delta t \left[ \mathbf{F}^n +  \boldsymbol{\nabla} H^n \right] \label{it:FSProject}
   \end{equation}

   The term $\hat{\mathbf{u}}^{n+1}=\mathbf{u}^{n} - \Delta t \mathbf{F}^n$ is known as intermediate velocity, and is a non matching divergence approximation to $\mathbf{u}^{n+1}$ (i.e. $\nabla \cdot \hat{\mathbf{u}}^{n+1} \neq \left( \nabla \cdot \mathbf{u}^{n+1} \right)^{th}$).
\end{enumerate}
%
Although in the original problem, Eqs.~\eqref{eq:LowMachMom}-\eqref{eq:LowMachDiv}, no boundary condition is required for $H(\mathbf{x},t)$, a consequence of the projection scheme is that boundary conditions are required on the Poisson equation of step Eq.~\eqref{it:FSPoisson}. For explicit methods and stationary \textit{solid} boundaries, the corresponding boundary condition is \textit{homogeneous} Neumann for $H^n$, $\partial H^n / \partial x_n =0$ in $\partial \Omega,\partial \Omega_1,...,\partial \Omega_{nbods}$ (i.e.~\cite{Perot:1993}).


The next component of the scheme involves approximating velocity boundary conditions for immersed solid boundaries. In FDS, velocity field components are staggered on Cartesian or cut faces. 
For complex geometry the projection scheme defined above is applied in the unstructured mesh composed by cut-cells and regular gas cells. Therefore, the Poisson equation is discretized in an unstructured grid and solved by matrix direct solvers, either by FDS mesh (\ct{ULMAT}) or globally (\ct{UGLMAT}). Currently the numerical scheme allows only one cut-cell/face Cartesian per cell/face, due to the current need to define a single background pressure per Cartesian cell. The blocking treatment of split cells is described in the Users guide. The projection step to obtain final velocities involves grid cell sizes in regular faces and centroid to centroid distance in cut-faces.

Cross velocities and tangent stresses and vorticity emanating from the immersed boundary are defined on Cartesian faces being intersected by the geometry surface and surrounding Cartesian edges. We use a wall modeled stress imposition method (STM) to this end described next.
 
\subsection{Stress Method: Tangential Velocities Estimated Through Cartesian Stresses}

In the discrete momentum equations, Eq.~\eqref{eq:LowMachMomEu}, the expression for $\mathbf{F^n}$ is:
%
\begin{equation}
   \mathbf{F}^n = -\left[ \mathbf{u} \times \boldsymbol{\omega} + \tilde{p} \nabla \left( \frac{1}{\rho} \right) 
   + \frac{1}{\rho} \left[ (\rho-\rho_0) \mathbf{g} + \nabla \cdot \tau \right]  \right]^n \label{eqn:momfn}
\end{equation} 
%
In this section we are interested in the divergence of deviatoric stress $ \nabla \cdot \tau$ term, its discretization and shear stress values at the boundary.

\subsubsection{Wall Modeled Shear Stress for Grid Aligned Solid}

% figure
\begin{figure}[h]
      \centering
       \includegraphics[trim = 70mm 60mm 70mm 50mm, clip,
       width=0.45\linewidth]{../../../fig/fds/Geometry_Figures/CutCellsWallAlignedSketch.png}
      \caption[Wall stress estimation]{Wall stress $\tau_w$ estimation using a wall function for solid boundary crossing through vertex $B$ and aligned with the fluid grid.}
        \label{Fig:Strss}
\end{figure}
%

Consider the grid aligned case of Fig.~\ref{Fig:Strss}. A velocity component $u_e^n$ is defined in a face adjacent to the solid boundary. 
The vertex point on the wall for this face is point $B$. A known shear stress $\tau_w$ is provided at this boundary point B. In two dimensions, the $x$ source component $\nabla \cdot \tau$ at point $e$ can be computed as:
%
\begin{equation}
   \left[ \frac{\partial \tau_{xx}}{\partial x} + \frac{\partial \tau_{xy}}{\partial y} \right]^n \simeq 
   \left[ \frac{\tau_{xx,E}-\tau_{xx,W}}{\Delta x} + \frac{\tau_{xy,N}-\tau_w}{\Delta y}   \right]^n
\end{equation}
%
Therefore, in an explicit time integration scheme, a boundary value of shear stress $\tau_w$ has to be estimated and used as source in the momentum equations. In FDS, this procedure is done using a wall model.
A simple wall model representation in this case is given by a function $u^+ = f(y^+)$, and the near wall scaling relationships:
%
\begin{equation}
  u_\tau = \sqrt{\frac{\tau_w}{\rho}}; \; \; \; y^+ = \frac{ y}{\delta_v}; \; \; \; \delta_v = \frac{\nu}{u_\tau}; \; \; \; u^+ = \frac{u}{u_\tau}
\end{equation}
%
Here $\tau_w$ is the stress at the wall point $B$; $\rho$ and $\nu$ are the density and molecular kinematic viscosity. The quantity $\delta_v$ is the near wall viscous length scale. The functional model used has been described in Sec.~\ref{info:velocity_bc}, and corresponds to a zero pressure gradient equilibrium turbulent boundary layer (\textit{log-law}). As this function is non-linear, $y^+$, $\tau_w$, $u_\tau$, and $\delta_v$ are computed to convergence through fixed point iteration. 

Referring to Fig.~\ref{Fig:Strss}, and considering the expression for the wall stress $\tau_w = \rho \nu \left( \frac{\partial u}{\partial y} \right)_{y=0}$ and a known velocity value $u_e^n$, initialize:
%
\begin{itemize}
   \item compute $\tau_w^0 = \rho \nu \frac{2 |u_e^n|}{\Delta y}$
   \item compute $u_\tau^0 = \sqrt{\tau_w^0/\rho}$
   \item compute $\delta_v^0 = \nu/u_\tau^0$
   \item compute $y^{+,0} = \frac{\Delta y}{2 \delta_v^0}$
\end{itemize}
%
If a direct numerical simulation (DNS) is performed, this is all that is needed to compute the wall stress (its definition). For wall modeled large eddy simulation (LES), for $k=1,2,..$ up to desired convergence do:
%
\begin{enumerate}
   \item compute $y^{+,k} = \frac{\Delta y}{2 \delta_v^{k-1}}$
   \item if $y^{+,k} < 12$ : compute $\tau_w^k = \rho \left( u_e^n / y^{+,k} \right)^2$; else use 
   $\tau_w^k =\rho \left( \frac{u_e^n}{\ln{(y^{+,k})}/\kappa + Br} \right)^2$
   \item compute $u_\tau^k = \sqrt{\tau_w^k/\rho}$
   \item compute $\delta_v^k = \nu/u_\tau^k$
\end{enumerate}
%
where  $\kappa \simeq 0.41$ and $Br \simeq 5.2$ are the log law constants. Once the stress at the wall $\tau_w$ is computed, it is used in the right hand side definition of the momentum equations to evolve the velocity flow variable $u_e^n \rightarrow u_e^{n+1}$. This velocity component
off the wall is integrated in time in the same manner as the rest of the discrete velocities, having to obey CFL and viscous constraints for stability.

\subsubsection{Non Grid Aligned Surfaces: Gradient Estimation in Cartesian Directions}

When the solid surface is not aligned with the fluid mesh, velocity gradients and Cartesian shear stress components are computed in the body surface and then are interpolated or extrapolated to surrounding Cartesian grid edges. In general, the approach taken is similar to what is presented in references~\cite{Lundquist:2010,MaLiu:2017}, but in our case, sampling points for the different Cartesian gradient components are defined in the Cartesian instead of wall normal directions.
In Fig.~\ref{Fig:STMgraphs}a, we define inside-boundary (IB) edges as edges contained completely in the solid (or solid surface) but adjacent to the boundary. Adjacency is defined as having one of the edge connected Cartesian faces being intersected by the solid boundary. We also define regular-cut face (RC) edges as edges completely defined in the gas phase, but adjacent to the boundary. Edges defined as IB or RC will receive vorticity and shear stress components extrapolated or interpolated from values estimated in the boundary through a wall model and gas phase values. The procedure might require an external gas phase edge (EP), also depicted in the figure.

% Figure
\begin{figure}[h]
      \centering
      \includegraphics[trim = 20mm 50mm 20mm 35mm, clip,
       width=0.95\linewidth]{../../../fig/fds/Geometry_Figures/CutCellsSTMgraphs.png}
      \put(-350,0){(a)}
      \put(-120,0){(b)}
      \caption[Velocity gradients and stress components]{Velocity gradients and stress components estimation. (a) For cut-face (i,j,k) velocity gradient $\partial u/\partial y$ and shear stress component $\mu \partial u/\partial y$ are estimated using a wall function for solid boundary crossing through vertex $B$, and transferred to Cartesian edges IB, RC. (b) Details of wall function used in obtaining Cartesian velocity derivative, shear stress component.}
        \label{Fig:STMgraphs}
\end{figure}
%

For each Cartesian face containing IB and RC edges, shear stress components are modified to account for the presence of the boundary.We assume the average location $B$ of the boundary is known in each direction for each Cartesian face containing cut-faces. In the following, all variables are assumed at time level $n$. In the case of $x$-face (i,j,k) in the figure, $\tau_{xy,IB}$ and $\tau_{xy,RC}$ are redefined as described in the following:
%
\begin{enumerate}
    \item For each RC edge compute the gas finite difference estimation of the cross velocity derivatives $\partial u_i/\partial x_j$ and $\partial u_j/\partial x_i$. In Fig.~\ref{Fig:STMgraphs}b, $\frac{\partial u}{\partial y} \simeq 2 \frac{u_E-u_F}{\Delta_y + \Delta_{RC}}$, where $u_E=u_{i,j+1,k}$ and $u_F$ is the cut-face velocity, and $\frac{\partial v}{\partial x} \simeq \frac{v_{i+1,j,k}-v_{i,j,k}}{\Delta_x}$. Note that the distance to the wall $\Delta_{RC}$ is being used to compute the finite difference in the cut-face and $y$ direction. Also interpolating the dynamic viscosity to the RC edge, estimate $\mu_{RC} \frac{\partial u}{\partial y}$ and $\mu_{RC} \frac{\partial v}{\partial x}$. \label{point1}
     \item Focusing in the face containing IB and RC edges, define point UB at a distance $\Delta y/2$ from the boundary B. Obtain velocity component $u_{UB}$ at location UB by linear interpolation of $u_E$ and $u_F$ in the $y$ direction.
     \item Using the wall function compute $\partial u/\partial y$ and $\mu_w \partial u/\partial y$ at boundary point B. Note that the input into the wall function routine is the velocity component $u_{UB}$ a distance $\Delta_y/2$ along the Cartesian direction $y$. We are using the wall model in a direction other than the wall normal direction. 
     \item The external edge EP is defined as: Take EP to be the RC edge if $\Delta_{IB} <= \Delta y/2$, use edge out in the Cartesian direction (i,j+1,k) otherwise. Compute $\partial u/\partial y$ and $\mu \partial u/\partial y$ at EP.
     \item Linearly extrapolate $\partial u/\partial y$ and $\mu \partial u/\partial y$ to IB edge from boundary B and EP values.
     \item If EP is the outer edge to RC, linearly interpolate $\partial u/\partial y$ and $\mu \partial u/\partial y$ to RC from EP and B values. To allow for $u_F$ feedback on the computation of this quantity we average the interpolated $\partial u/\partial y$ and $\mu \partial u/\partial y$ with the value computed in point~\ref{point1}.
\end{enumerate}
%

Finally, all RC and IB edges are assigned a vorticity and shear stress using the computed gradients $\partial u_i/\partial x_j$ and components $\mu \partial u_i/\partial x_j$. The computed vorticity and stresses are used in the construction of the force vector $\mathbf{F}^n$ of equation~\eqref{eqn:momfn} and used in the projection equations~\eqref{it:FSPoisson}-\eqref{it:FSProject}.


In similar manner as in scalar transport, small cut-faces and their associated volumes will constraint the momentum explicit time integration steps for stability. A similar linking procedure as described for scalars is employed to link small cut-faces to surrounding faces and to define a single average value of $\overline{\mathbf{F}}^n$ for the linked set. Pressure gradient values on these faces are not linked, as required for correct cell by cell enforcement of the discrete divergence constraint and mass conservation.

\section{Energy Equation, Thermodynamic Divergence}

\subsection{Divergence Constraint}

Starting from the sensible enthalpy evolution equation, the divergence of the velocity field  imposed in eq.~\eqref{eq:LowMachDiv} can be factored as:
%
\begin{eqnarray}
    ( \nabla \cdot \mathbf{u} )^{th} &=&
    \left[ \frac{1}{\rho c_p T} - \frac{1}{\bar{p}} \right]
    \frac{\partial \bar{p}}{\partial t} + \frac{w \rho_0 g_z}{\rho c_p T} \nonumber \\
    &+& \frac{1}{\rho c_p T} \left[ \dot{q}''' - \nabla \cdot \dot{\mathbf{q}}'' - \mathbf{u} \cdot \nabla (\rho h_s) \right] \nonumber \\
    &+& \frac{1}{\rho} \sum_\alpha \left( \frac{\overline{W}}{W_\alpha} - \frac{h_{s,\alpha}}{c_p T} \right) \left[ \dot{m}_\alpha''' - \nabla \cdot \mathbf{J}_\alpha - \mathbf{u} \cdot \nabla (\rho Y_\alpha) \right] \label{eq:divth}
\end{eqnarray}
%
we call this divergence expression the thermodynamic divergence $( \nabla \cdot \mathbf{u} )^{th}$. The projection scheme for velocities enforces this final divergence on the discrete velocity field in each cell of the spatial discretization.
Next, we look at the finite volume discretization of the non-conservative advection terms involving the velocity $\mathbf{u}$, in the right hand side of Eq.~\eqref{eq:divth}.

We will use the finite volume version of the previous equation to obtain the volume integrated target thermodynamic divergence on each cut cell volume, and to recompute the thermodynamic divergence on regular gas phase cells adjacent to cut-cells.


\subsection{Finite-Volume Discretization of Scalar Advection Terms}

Consider the the equality
%
\begin{equation}
    \nabla \cdot \mathbf{\phi u} = \mathbf{u} \cdot \nabla (\phi) + \phi \nabla \cdot \mathbf{u} \label{eq:advforms}
\end{equation}
%
Integrating over the volume of cell $ii$
%
\begin{equation}
    \int_{\Omega_{ii}} {\mathbf{u} \cdot \nabla(\phi)} d\Omega =
    \int_{\Omega_{ii}} \nabla \cdot ({\phi} \mathbf{u}) d\Omega -
    \int_{\Omega_{ii}} \phi \nabla \cdot \mathbf{u} d\Omega \label{eq:flxlim1} \\
\end{equation}
Assuming a flux limited interpolation of the scalar $\phi$ at the cell boundaries, the discrete counterpart is
%
\begin{equation}
    \overline{\mathbf{u} \cdot \nabla(\phi)} \; V_{ii} =
    \sum_{k=1}^{nf_c} (\overline{\phi} \mathbf{u})_{ii,k} \cdot \hat{\mathbf{n}}_{ii,k} \: A_k -
    [\phi]_{ii} \sum_{k=1}^{nf_c} \mathbf{u}_{ii,k} \cdot \hat{\mathbf{n}}_{ii,k} \: A_k \label{eq:flxlim2}
\end{equation}
%
Here the overline states that the cell-centered scalar quantity $\phi \in \{\rho h_s, \rho Y_\alpha\}$ has been interpolated to the cell faces using a flux limited interpolation scheme. The overline in $\overline{\mathbf{u} \cdot \nabla(\phi)}$, states that this term is computed consistent with the flux limited interpolation adopted. Also, $V_{ii}$ is the volume of cell $ii$, and $A_k$ is the area of a face $k=1,...,nf_c$ of the given cell. This last equation is the finite volume counterpart of Eq.~\eqref{eqn_flux_decomposition}. %(B.12).

\subsection{Discrete Thermodynamic Divergence Expression}

For cut-cell $ii$, the corresponding discrete volume integrated expression is:
%
\begin{eqnarray}
    ( \nabla \cdot \mathbf{u} )_{ii}^{th} \; V_{ii} &=&
    \left[ \frac{1}{(\rho c_p T)_{ii}} - \frac{1}{\bar{p}_{ii}} \right]
    \frac{\partial \bar{p}_{ii}}{\partial t} \; V_{ii} +
    \frac{w_{ii} \: \rho_0 g_z}{(\rho c_p T)_{ii}} V_{ii} \nonumber \\
    &+& \frac{1}{(\rho c_p T)_{ii}} \left[ \dot{q}''' V_{ii} - \nabla \cdot \dot{\mathbf{q}}''_r V_{ii} -
    \sum_{k=1}^{nf_c} \dot{\mathbf{q}}''_{ii,k} \cdot \hat{\mathbf{n}}_{ii,k} \: A_k
    - \overline{\mathbf{u} \cdot \nabla (\rho h_s)} V_{ii} \right] \nonumber \\
    &+& \frac{1}{\rho_{ii}} \sum_\alpha \left( \frac{\overline{W}}{W_\alpha} - \frac{h_{s,\alpha}}{c_p T} \right)_{ii} \left[ \dot{m}_\alpha''' V_{ii} -
    \sum_{k=1}^{nf_c} \mathbf{J}_{\alpha,ii,k} \cdot \hat{\mathbf{n}}_{ii,k} \: A_k
    - \overline{\mathbf{u} \cdot \nabla (\rho Y_\alpha)} V_{ii} \right] \label{eq:divth2}
\end{eqnarray}
%
where the overline terms refer to flux limited interpolation of corresponding scalars, terms defined with subscript $ii$ refer to cell defined quantities, and the scalar diffusive flux $\mathbf{J}_\alpha=- \rho D_\alpha \nabla Y_\alpha$. The divergence of radiative heat flux $\nabla \cdot \dot{\mathbf{q}}''_r$ is provided by the solution of radiative transport equation, Sec.~\ref{radnumericalmethodsection}.  Also, the vertical velocity $w_{ii}$ \textit{has been interpolated} to the cut-cell centroid. The flux vector $\dot{\mathbf{q}}''_{ii,k}$ refers to diffusive and conductive heat fluxes described next.

\subsection{Computing Diffusive Heat Flux in Finite-Volume Form}

Given the species diffusion model adopted, the diffusive heat flux vector field is $\dot{\mathbf{q}}_{d,\alpha}''=-h_s \rho D_\alpha \nabla(Y_\alpha)$. The integral of its divergence over a cell $ii$ is:
\begin{equation}
    \int_{\Omega_{ii}} {\nabla \cdot \left(-h_{s,\alpha} \rho D_\alpha \nabla(Y_\alpha) \right)} d\Omega = \sum_{k=1}^{nf_c} \left(-h_{s,\alpha} \rho D_\alpha \nabla(Y_\alpha) \right)_{ii,k} \cdot \hat{\mathbf{n}}_{ii,k}  \: A_k
\end{equation}
where $\left(-h_s \rho D_\alpha \nabla(Y_\alpha) \right)_{ii,k}$ is the mean heat flux defined on boundary surface $k$ of cell $ii$.


\subsection{Computing Heat Conduction Flux in Finite-Volume Form}


The heat conduction flux vector of the mixture is defined as $\dot{\mathbf{q}}_{kT}''= -k \nabla(T)$, where $k$ is the mixture conductivity. The integral of its divergence over a cell $ii$ is computed similarly as in the previous section:
\begin{equation}
    \int_{\Omega_{ii}} {\nabla \cdot \left(-k \nabla(T)\right)} d\Omega = \sum_{k=1}^{nf_c} \left(-k \nabla(T) \right)_{ii,k} \cdot \hat{\mathbf{n}}_{ii,k}  \: A_k
\end{equation}
where $\left(-k \nabla(T) \right)_{ii,k}$ is the mean heat flux defined on boundary surface $k$ of cell $ii$.


